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ABSTRACT 


Various approaches to the frequency-sampling design of two-dimensional FIR 
filters are analyzed. The IDF T approach requiring uniform sampling on a Cartesian 
grid is first described. A method which allows arbitrary placement of frequency 
samples but which does not satisfy the Haar condition is presented. Finally, a 
novel, computationally efficient method which allows nonuniform sampling and 
which always provides a unique design solution is presented. The new approach 
is compared with the other methods in terms of design flexibility, computational 
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I. INTRODUCTION 


Finite impulse response (FIR) filters have seen widespread use in the field of 
two-dimensional digital signal processing. Reasons for this include the inherent 
stability of FIR filters, the ability to achieve a linear (or zero) phase characteristic 
in the frequency response, and relative ease of design. Traditional approaches to 
the design of such filters include the window method, in which a smoothing window 
is applied to the Fourier series coefficients of the ideal frequency response; the fre- 
quency transformation method, in which a one-dimensional (1-D) prototype filter 
response is mapped into a function of two frequency variables using an appropri- 
ate transformation function; and the frequency—sampling method, in which desired 
frequency response values are specified at certain sample points in the frequency 
domain. Of the three methods described above, the frequency-sampling method is 
probably the least utilized and the least understood, and hence warrants further 
investigation. This thesis deals with a comparative analysis of various methods of 
designing two-dimensional (2-D) FIR filters using frequency-sampling techniques. 
In the analysis, some new methods for frequency-sampling design are presented. 

The 2—D frequency-sampling design method is, in essence, an interpolation 
problem in two variables where the two variables are the two frequencies of interest, 
denoted as w, and wz. The frequency response is specified at selected sample 
points in the (w,,w2) plane, and is expressed in terms of linearly independent 
basis functions. The solution to this interpolation problem involves finding the 
coefficients associated with the respective basis functions in order to meet the 
interpolating conditions. These coefficients are, in general, closely related to the 


impulse response. The various approaches to the frequency-sampling design 


problem involve differing constraints placed on sample location in the (w, ,w,) plane 
and differing sets of basis functions. Because the resulting frequency response takes 
on exactly the prescribed value at each frequency sample location, the response is 
somewhat predictable at the onset of the design process. A disadvantage of 2-— 
D frequency-sampling design, however, is that the existence and uniqueness of a 
solution is dependent upon the location of selected frequency samples and the basis 
functions specified; in other words, degenerate cases may arise in which no solution 
exists to the bivariate interpolation problem. 

FIR filter design techniques denoted as “frequency—-sampling” traditionally in- 
volve sampling a desired frequency response at equally-spaced intervals and then 
solving for the impulse response coefficients using an inverse discrete Fourier trans- 
form (IDFT) [Ref. 1]. For the 2—-D case, the process is similar, with the desired 
frequency response sampled at equally spaced points on a Cartesian grid in the 
(w,,w2) plane [Refs. 2, 3]. This approach will be referred to as the “uniform 
sampling” method in this thesis. There are two significant benefits to such an ap- 
proach in two dimensions. First, with frequency samples placed at the vertices of 
an N x N Cartesian grid, the existence of a unique filter satisfying the interpolating 
conditions and possessing an impulse response with an N x WN region of support is 
guaranteed. Second, use of the IDFT, or alternatively, the fast Fourier transform 
(FFT) algorithm, assures a computationally efficient design process. The princi- 
pal disadvantage of the approach is that there is no flexibility in the placement of 
frequency samples. The approach is presented in more detail in Chapter II. 

A different approach to frequency-sampling design of 2-D FIR filters involves 
allowing frequency samples to occur anywhere in the (w,,w,) plane and fixing 
the region of support for the filter impulse response. This approach is presented 


in Chapter III and will be denoted as the “arbitrary sampling” technique. The 


method is quite flexible since samples can be arranged to enhance certain desired 
features in the resulting frequency response, such as flat passband, sharp roll-off, 
equi-ripple behavior, etc. The design, however, involves the solution of a large 
system of linear equations for a filter of any practical order and hence it is quite 
computationally intensive. Additionally, with no constraints on the location of 
samples, degenerate cases may arise in which no solution exists which meets all 
interpolating conditions for the specified impulse response region of support. 

A novel 2—-D frequency—sampling design technique is presented in Chapter IV, 
and is referred to as the “nonuniform sampling” method. This method places 
certain constraints on frequency sample location but still possesses much greater 
flexibility than the uniform sampling approach. Rather than solving one large 
system of equations, this technique involves the solution of several smaller systems 
of linear equations. Given the constraints on the sampling geometry, existence and 
uniqueness of a solution is guaranteed. Additionally, the design method is less 
computationally intensive than the arbitrary sampling case but more intensive, in 
general, than the uniform sampling case. 

In comparing and contrasting these various approaches to 2—D frequency- 
sampling design, several interesting points arise. First, frequency-sampling design 
techniques involve some sort of decision-making process in determining where to 
place the samples in order to best match the desired frequency response. Methods 
which are more constrained, such as the uniform sampling approach, involve a lesser 
degree of decision-making and hence can be considered easier to use from that 
perspective. Second, by reducing the number of constraints on the placement of 
frequency samples, flexibility in matching the desired frequency response is gained. 
Finally, and perhaps most significantly, computational efficiency in determining 


filter parameters can be directly related to the number of constraints on sample 


location. Of the techniques discussed here, those which are more constrained tend 
to be more computationally efficient. Thus there are tradeoffs to be considered 
when selecting one of these techniques. In this thesis, a detailed analysis of these 
concepts is performed and the quality of these frequency-sampling methods in 


matching desired filter characteristics is investigated through design examples. 
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II. DESIGNS WITH UNIFORM SAMPLING ON A CARTESIAN GRID 


The approach to frequency-sampling design of 2—-D FIR filters discussed in this 
chapter is the one referred to as “frequency-sampling” in various texts [Refs. 1, 
2|. It will be referred to the “uniform sampling” method through the remainder of 
this thesis. This method has significant computational advantages through use of 
an inverse discrete Fourier transform (IDFT) approach. Additionally, this method 
will never result in any degenerate cases, i.e. situations in which a filter satisfying 
each of the interpolating conditions cannot be found. However, due to constraints 


on the location of samples, the method is rather inflexible. 


A. APPROACH 
The transfer function of a two-dimensional FIR filter with impulse response 


support in the first quadrant can be defined as 


1-~1N3-1 


gee. | + S- AL (use |e ey (2.1) 


mi=O0 n4=0 
The frequency response of a such a filter is obtained by evaluating the transfer 
function on the unit bicircle through the substitutions z, = e?“! and z, = e?”?, 


yielding 


N;-1N3,-1 


A (wy, we) = yh ilgite |e 2 ee (2.2) 


Np 0 Ns 0 


The frequency-sampling filter design problem involves determination of a sequence 
h(n, ,n2) which will result in a filter matching a desired frequency response H,(w,,w2) 
at specified sample points in the (w,,w 2) plane. 

The initial step in the design process is to choose the filter order, and hence 


the parameters N, and N2. Often it is desirable to specify N, = Nz = N to 
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ensure a result with equal dependence on the frequencies w, and w,. The desired 
frequency response H,(w,,w 2) is sampled at the vertices of a uniformly spaced 
N, x N, Cartesian grid on the region {0 < w, < 27; O< we. < 27} of the (w,, we) 
plane. The resulting frequency samples can be expressed in 2—D sequence form as 


~ 


H(k,, kz) = Hy(w, Wa) | tees deka 
ky Sei leecee yal and cen). cea le (2.3) 


The resulting sampling geometry is illustrated in Figure 2.1. 


U4 


Ny, 








0 2n 4n 2n(N,-1) 








Figure 2.1. Arrangement of Samples for Uniform Sampling Method. 


Because an IDFT approach is to be used in computing h(n,,n,), and since the 
IDFT is defined for first-quadrant sequences only, the desired frequency response 


should possess a linear phase characteristic such that 
Hi(@,, Oo) = dg ois ee een eee am (2.4) 
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where 


Mi= {eye Nod (28) 
and 
— a8 a /2 Ny odd ae 
Additionally, H,(w,,w2) must possess the symmetry conditions 
\H(ki,k2) =|H(k,,N2 —ke)| 
Ca) 


= | (N, a k,,k2)| 


in order to ensure a purely real impulse response. 
The next step is to compute h(n,,n,) by taking the IDFT of the sampled 


desired frequency response: 


lige) = _ » ke) 


wo 2 ) H(k, k,) pee WT? 28) 
b ] 
k,=0 kz=0 
where 
—j7it 7 a 
Wee i Cand =) pee 2 "> . 


The resulting sequence h(n,,n-2) is purely real and represents the impulse response 
of the designed filter. Such a filter possesses a linear phase characteristic. A zero- 
phase characteristic can be obtained by shifting the impulse response so that it is 


centered about the origin such that 


h'(m,,m2) = h(m, + M,,m, + M2), (2.9) 


where h'(m,,mz2) represents the zero-phase impulse response. 
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To see the relationship between the resulting filter frequency response and the 


specified sample sequence H(k,,k2), Eq. (2.8) is substituted into Eq. (2.2): 





x le N,-1N,-1 


H(w1,W2) = y yy NN ss SH (ky, , ke eee 
2 


my= On >=0 k,;=0 ko= 0 
Ni,-1N},-1 N,-1 N 3-1 


= ra > SS H (ky, kz) Ss. yee 223) >. (Wy. 67a 


k,=0 k2=0 nm ,=0 maz=0 





-> bs (1-—e7 aN Vl er ae H(k1, kz) 
Te Sa N,N, (L—W,,*e-##:)(1-—Wy*te-fea) 
(2.10) 


The above expression can be rewritten in the form 


N,-1N3-1 


Ay 
H(w,,we) = Sk CAL, i Wy et w2 am wy, ka)» (2.11) 
1 2 


EOE 0 
where 


(1 = ee al a emia 2)) 


O(w, , We ) ~ N,N, (1 — e-2#)(1 — e242) 
: p27 991 (4) eee = Wy) er? ™ea(t- ee) 


which simplifies to 


sin(w, N, /2) sin(w. N, /2) 
N, Nz sin(w, /2) sin(w, /2) (2312) 


vero — jw3(—4+—) Apt as Wy) eitha(l- wy) 


®(w,, We) = 


for the case where N, and N, are both odd integers. This analysis is a 2—D 
extension of operations performed for the one-dimensional case [Ref. 1: pp. 97- 


98]. Since 


ee ee (k,, kz) _ (0, 0) 
(Fea seks) = 4 4 ky = 1,2; .<2,.N) = iieand he oe 


the frequency response is matched exactly at the original sample points. That is, 


H(w1,2)|,,_ 2281 4, — 282 = H(k,, ka), (2.13) 


ce Oey = le and k= 0,1,...,N,—1. 


Thus, the design process entails a bivariate interpolation in the variables w, and 
W2, and, due to use of the IDFT algorithm, the interpolating functions are of the 


Homma of Kg. (2.12). 


In using this filter design method, the choice of the desired frequency response 
H,(W,,w2) has significant effects on the quality of the resulting filter. For instance, 
if an “ideal” lowpass filter characteristic with the circular symmetry shown in 
Figure 2.2 is used for H4(w,,w.), the passband and stopband of the resulting filter 


frequency response exhibit substantial ripple. 
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Figure 2.2. Ideal Lowpass Filter Frequency Response. 
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Figure 2.3 illustrates contour and perspective plots of the frequency response 
arising from such a design. In this case, a 15 x 15 point filter was designed by 
uniformly sampling the desired frequency response of Figure 2.2. 

The frequency response of the resulting filter can be improved considerably if 
a transition region is specified for Hy(w,,w2). Increasing the width of such a tran- 
sition region reduces passband and stopband ripple but does so at the expense of a 
less sharp roll-off characteristic in the resulting filter response. The values assigned 
to H (ky ,k2) for points (k,,k2) located in this transition region significantly influ- 
ence the behavior of the resulting filter. The selection of these values in order to 
satisfy some specific optimization criterion for the frequency response, such as best 
equi-ripple approximation to an ideal response, is not trivial. Since H(w,,w2) isa 
linear function of the values H(k, ,k2), as seen in Eq. (2.10), linear optimization 
techniques can be employed in selecting the transition region values of H(k; Ke). 
Such techniques for the one-dimensional case are rather involved [Ref. 4]; for the 
2—D case, the problem is much more computationally intensive [Ref. 2], and hence, 
has not been addressed. Specifying H,(w,,w2) to have a smooth, linear transition 
between the passband and stopband regions is the simplest thing to do and gen- 
erally works well. The magnitude characteristic of a desired frequency response 
with a linear, circularly symmetric transition region is shown in Figure 2.4. A fil- 
ter can then be designed by sampling this response on the uniform Cartesian grid 
specified in Eq. (2.3). Contour and perspective plots of a 15 x 15 filter designed 
in this manner are shown in Figure 2.5. In comparing this result to that of Figure 
2.3, it is evident that inclusion of a transition region greatly reduces passband and 


stopband ripple. 
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Figure 2.3. Frequency Response Arising From Sampling Ideal Lowpass 


Characteristic With No Transition Band. 
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H (w1,w2) 





Figure 2.4. Desired Lowpass Frequency Response Characteristic 
With Linear Transition Band. 


B. EXISTENCE AND UNIQUENESS 


A significant property of the uniform sampling design is that a unique filter 
satisfying each of the interpolating conditions specified in Eq. (2.3) is guaranteed 
to exist. In general, existence and uniqueness of such an interpolation result in two 
dimensions is not guaranteed; this can be a significant theoretical and practical 
problem. The point is discussed in further detail in Chapter III. A proof of existence 
of a unique h(n,,n,) using the frequency-sampling method with uniformly spaced 


samples follows. 
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Figure 2.5. Frequency Response Arising From Sampling Desired Lowpass 
Characteristic With Linear Slope Transition Band. 
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Equation (2.8) is used to determine h(n,,n2). This can be rewritten as 


ae Ne 
il 1 - n = ni 
Gori. ee | = aE ye Fe ys H(k,,k, Ae nia Wee re (2.14) 
ky, = 0 k,=0 
Now define 
ce 
G(k,,n2) = re ae ie (Zale) 
k3,=0 


Substituting Eq. (2.15) into Eq. (2.14), the following expression for the impulse 


response results: 
Nae 


1 
h(n, ,n2) => N, xe (k,, 72) We AED (2.16) 


Thus the 2-D IDFT of H (ky ,k,) can be determined in a two-step process involving 
one-dimensional transforms known as row-column decomposition |Ref. 5: pp. 75- 
76]. 

The first step in the row-column decomposition process is to perform the 
solution of a total of N, systems of N. equations of the form of Eq. (2.15) using a 


1-D IDFT algorithm. In matrix form, this can be written as 


G(k; ,0) ee Vie Ra Wr, H(k, 0) 
G(k,1) eV ce We ase H (ky ,1) 
: Na |: ee | | 
G(k,,N2 — 1) WoW Oe en H(k,,N2 — 1) 
(2.17) 


which can be expressed in vector-matrix notation as 


1 fe 
G = —- Wil. (2.18) 


2 


It can be shown [Ref. 6: pp. 44-45] that a scaled version of W,, namely W2/V Nz, 
is unitary and hence non-singular. Therefore a unique solution for G exists for each 


k;, p—tOet,...,N, — 1. 
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The second step in the computation involves the solution of N, systems of N, 


equations of the form of Eq. (2.16). In matrix form, this becomes 


h(0, no) Wir DAs Sagat: AAS G(0, nz) 
him) | a [wee wet we] Ginn) 
3 eel alen: ae 3 | 
h(N, — 1,72) 1 ella | aah a Ware G(N, — 1,72) 
(2.19) 


which in vector-matrix form is 


ae 


h 
N, 


WG. (2.20) 


Since G was uniquely determined for each value of k,, and since W, is also 
non-singular, a unique solution for h exists for each n;, jg = 0,1,...,N. —1. 
Therefore, a unique solution for h(n,,n) is obtained from the sequence H (ky nk) 3 
which represents the samples of a desired frequency response at the vertices of a 


uniform Cartesian grid in the (w,,w,) plane. 


C. COMPUTATIONS 

While there is no flexibility in the choice of sample point locations using this 
IDFT approach, the method does offer some significant computational advantages. 
One criterion often used to measure computational efficiency in solving a problem 
is to determine the number of complex operations involved. While there exist many 
methods of computing the 2-D IDFT, the use of row-column decomposition and 
the IDFT directly (vice the FFT algorithm) is a reasonably efficient approach for 
the design of filters of practical order and is used in the following analysis. 

The initial step in the row-column decomposition approach is to solve a total 
of N, systems of the form of Eq. (2.17). Solution of each system via the IDFT 


involves a total of N? multiplications and additions. Therefore, a total of N, N3 
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floating point operations are required in this step, where a floating point operation 
is considered one multiplication and one addition. 

The second part of the row-column decomposition approach involves the so- 
lution of N, systems of the form of Eq. (2.19). Solution of each system via the 
IDFT requires N? operations. Therefore, N? N2 floating operations are involved 
in the second step. 

The total number of floating point operations required in the solution for the 
impulse response coefficients of an N, x N, frequency-sampling FIR filter design 
using uniform sampling on a Cartesian grid and a row-column decomposition, IDF T 
approach, is 


C inion = N, N: ais N? No. (2.21) 
For the case in which N, = N, = N, this reduces to 
Cinitoemn = ZINE (2e22) 


In summary, the uniform sampling, IDF T approach to 2—D frequency sampling 
filter design has the benefit of computational efficiency, yet lacks flexibility in the 
choice of frequency samples. The next two chapters will present methods which 


achieve greater flexibility by allowing nonuniform spacing of samples. 
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Il. DESIGNS WITH ARBITRARILY PLACED SAMPLES 


In Chapter II a frequency-sampling 2—D FIR filter design method was pre- 
sented in which samples were placed on a uniform Cartesian grid in the (w,,w2) 
plane. Due to the constraints on sample placement, however, the method was not 
particularly flexible. If the desired frequency response H,(w,,W-2) was fixed, no free 
parameters were available in optimizing the frequency response of the resulting fil- 
ter. In this chapter, an alternate frequency-sampling design technique for 2-D FIR 
filters is presented in which samples are no longer constrained to a uniform Carte- 
sian grid but can be placed anywhere in the (w,,w2) plane. This approach will 
be referred to as the “arbitrary-sampling” method. Allowing the samples to occur 
anywhere in the (w,,w,) plane introduces additional free parameters; hence, linear 
optimization techniques can be employed without varying H,(w,,w2), the desired 
frequency response to be sampled. Since the locations of the frequency samples are 
much less constrained, the method provides greater potential for achieving a desir- 
able frequency response characteristic in the resulting filter. Since the placement 
of samples on a uniform Cartesian grid is a special case of the arbitrary-sampling 
method, the arbitrary-sampling method can achieve results at least as good as 
the uniform-sampling method for a fixed impulse response region of support and 
a fixed number of samples. However, the arbitrary-sampling method possesses a 
few significant drawbacks, such as computational complexity and the possibility of 
degenerate cases in which the interpolating conditions cannot be met. These must 


be considered in the selection of a method for filter design. 
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A. APPROACH 

The arbitrary-sampling method can be used to produce 2-D FIR filters with 
a variety of impulse response regions of support. This region of support is one of 
the parameters to be specified at the outset of the design process. For illustration 
purposes, a filter with an N, x Nz. rectangularly-shaped impulse response region 
of support centered on the origin of the (n,,n,) plane is specified, where N, and 


N, are odd integers. The transfer function of such a filter can be expressed as 


My M; 
Eee | = SS MS Digiiewice | Zee (3.1) 


ny=—-M,n3,=-M3 


where h(n,,n2) is the impulse response to be determined, M, = (N, — 1)/2 and 
M, = (N, —1)/2. The parameters M, and Mp, will be used often in the remainder 
of this analysis, so the above relationships are worth noting. Evaluation of Eq. 


(3.1) on the unit bicircle yields the following expression for the frequency response: 


M, M; 
cme — S- SS ientg ye" ane tam (s22)} 


nmy=~-M, nyz3=—-M, 


To ensure a zero-phase characteristic in the resulting frequency response, 


h(n,,n2) is specified to possess quadrant symmetry such that 
h(ni,n2) = h(—n,,n,) = h(n, —nz2). (3.3} 


Application of Eq. (3.3) to Eq. (3.2) produces the following expression: 


My, 

H(w,,w2) =h(0,0) + > 2h(n,,0) coswyn,+ 
my=l 

My, M, M3 


S- 2h(0, nz) cos W2n2 + » << 4h(n,,n2) cosw,n, cos W2N2. 


n5=1 my=ln4,=1 


(3.4) 
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This can be rewritten as 


M, M32 


H(w,,w2) = ‘A Sy a(n, ,N2) COS W1N, COSW2N2, (3.5) 


na —v i> 0 


where 
h(ny,n2) nm, =0,n2 =0 


2h(n,,n-) Nj # 0,n- = 0 


= 3.6 
a(n ,72) 2h leet — Os, 0 ee) 
AT Guede ie 70,72, 7- © 

The basis functions {w~,;(w,,w2)} are defined such that 

Wb; (Wi, W2) = Cosw, nN, COSW2N2, (3.7) 
where . 
1=n,(M,4+1)+n, +1; ee Oe IVE: 

(3.8) 


No = OP ee iis. 
By arranging the 2—D sequence a(n,,n-2) into the one-dimensional sequence a(t), 


the frequency response can then be expressed as 


(co ea) = Ss a(t); (wi ,W2), (3.9) 


R 
ell 


1 
where 


R=(M, +1)(M, +1). (3.10) 


From Eq. (3.9), the elements of the sequence a(z) represent a total of R free 
parameters to be determined. Therefore, the frequency response is to be specified 
at a total of R distinct sample points in the (w,,w,) plane. 

The initial step in the design process is to specify H(w,,w,) at R distinct 
sample points in the region K = {(w,,w2.):O0<w, <7;0<w, < 7}. An example 
of such a sampling geometry is illustrated in Figure 3.1. 

No constraints other than the distinctness requirement are placed on the loca- 


tion of samples within this region. The form of Eq. (3.5) implies that the frequency 
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Figure 3.1. Possible Arrangement of Samples for Arbitrary Sampling Design. 


response possesses quadrantal symmetry and periodicity of 27 in both frequencies. 


Hence, H(w,,w.) is defined for all (w,,w2). 
The next step in the design process is to solve for the sequence a(n, ,n2). Appli- 
cation of Eq. (3.9) at each specified sample location {(w,,w2); ; «= 1,2,...,R} 


results in the following system of linear equations: 


H(w,,we), Vv, CARO the (Wi We), eatin (w,,W2), a(1) 

H (wy ,W2)2 a vy, (w,,We)2 2 (W, , We )2 eR (W, ,W2)2 a(2) 

Hien Y, (reais Wp. (Wi ,Wo)p semis ir (W,,W2)r a(R) 
(sam 
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In vector-matrix notation, Eq. (3.11) becomes 
lef == 181, (3212) 


in which B is of dimension R x R. The sequence a(n, ,n-2) results from solving the 


relation 


a= > He (8.13) 


The final step in the process is to compute the filter impulse response h(n, ,n2) 
using the relationship of Eq. (3.6). The resulting impulse response has an N, x N, 
region of support. Equation (3.4) suggests a transfer function of the form 


My, M3 
Mere 2 (0,0) Sy t(n,,0)(2)° 4+ 2," Je * A(O,n2) (227 + 2, "") 


m,=1 m4o=—1 


My, M 4 


+ SO YE Afr may(et + 2p" (ep? +27"); 

I (3.14) 
which can be used in implementing a 2-D filter structure with a minimal number 
of gains. 

This same approach can also be employed if other than rectangular impulse 
response regions of support are specified. Again, a system of linear equations of 
the form of Eq. (3.5) is solved, but the limits on the summations are altered to 
reflect the desired region of support for h(n,,nz). 

To summarize the arbitrary-sampling design process, the following steps are 


performed: 


(1) Fix the impulse response region of support, and hence, the filter order. 


(2) Specify the desired frequency response at selected sample points on the 
region K = {(w,,w2):0<w, <7;0<w. <a}. The number of specified 
samples must equal the number of independent impulse response values 
(i.e., those which are not fixed by the symmetry conditions). The 


uf 


locations of samples within the region is unconstrained except for 


the condition that each sample location be distinct. 


(3) Solve a system of linear equations to obtain the resulting impulse 
response h(n, ,n2). 


B. EXISTENCE AND UNIQUENESS 

There is a fundamental deficiency in the arbitrary-sampling approach to the 2- 
D FIR frequency-sampling filter design problem. For specified frequency response 
samples placed arbitrarily in region K, there is no guarantee that a unique filter 
possessing a number of independent impulse response samples equal to the number 
of interpolating conditions, and which satisfies each of those interpolating condi- 
tions, exists. This is a problem inherent to interpolation in two dimensions which 
does not arise in the one-dimensional case. 

For a unique solution to the system of equations indicated by Eq. (3.12) to 
exist, the basis functions {,;(w,,w2)} must satisfy the Haar condition. A set of 
vectors in n-space satisfies the Haar condition if every set of n of them is linearly 
independent |Ref. 7: p. 45]. Therefore, to satisfy the Haar condition, any set of 
R characteristic vectors defined on RF distinct points in region K of the (w,,w2) 
plane must be linearly independent. For the arbitrary-sampling method developed 
in this chapter, a characteristic vector {W((w,,w2);)} associated with a particular 


sample location (w,,w2); can be defined as 


W((w1,W2)i) = [Pi ((w1,W2)i) Po((wi,We)i) --- Dr ((wi,2):)]° - (3.15) 


It is known that for a region K' of the two-dimensional frequency plane, there are no 
nontrivial sets of two-dimensional basis functions which satisfy the Haar condition 
[Ref. 8: p. 494]. Therefore, there is no guarantee that the sets of characteristic 


vectors formed by evaluating Eq. (3.15) at arbitrarily selected sample points in 
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region K of the (w,,w2) plane are linearly independent. Since the square matrix 


B in Eq. (3.12) can be expressed in terms of the characteristic vectors as 


Mes a OMe 6 U((w.e)n) | (3.16) 


and since the characteristic vectors may be linearly dependent, B may be singular 
and a solution to Eq. (3.13) may not exist. Such a case is termed a degeneracy, 
and illustrates a basic problem associated with the arbitrary-sampling method. As 
can be inferred from the form of the basis functions, the occurrence of a degeneracy 
is dependent solely upon the location of samples; the specified frequency response 
values at the sample points do not affect this. 
A simple example illustrates how a degeneracy can arise. The frequency re- 
sponse of a filter is specified at four sample points in the (w,,w.) plane: 
eer (Os0:67) =" 1 2. H(0.47,0) =1 
So (OAnea, = 0 4. H(x,0.67) =0 
These sample locations are shown in Figure 3.2. 
The parameters M, and M, are both set equal to 1 so that the number of 
samples in the sequence a(n,,n,.) equals the number of frequency samples. Sub- 


stitution into Eq. (3.11) yields 


1 1 cos0.67 1 cos 0.67 a(0, 0) 

IL Joc | at a cos0.4m  cos0.4n a(0, 1) (3.17) 
0 1 —1 cos0.4% —cos0.4m |} | a(1,0) | © 

0 1 cos0.67 —1 — cos 0.67 Aol 


Some calculations reveal that the square matrix above has a determinant of zero; 
therefore, it is not invertible and no unique solution for a(n,,n2) exists. 

To determine whether degenerate cases present a serious shortcoming of this 
method, it is of interest to determine how often degeneracies are likely to arise. To 


investigate this, an APL computer program was developed to count the number of 
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Figure 3.2. Arrangement of Samples for Design Example. 


degenerate cases arising from the arbitrary selection of sample points {(w,,w2); } 
on region K of the (w,,w2) plane. A total of R samples were arbitrarily placed 
at distinctly different vertices of a grid of specified density on region K for each 
case. The placement of samples was performed in a pseudo-random fashion using 
the APL “roll” function. The density of the grid, in effect, is determined by the 
number of decimal places carried in the frequency values. In the analysis, two 
parameters were varied: the dimensions of the specified filter impulse response, 
and the grid density. For each situation, one hundred separate sets of arbitrary 
samples were investigated. The program then determined the number of cases in 
which the sampling geometry resulted in a degeneracy. Due to the large number of 
computations involved, analysis was limited to filters with square impulse response 


regions of support only. The results are listed in Table 3-1. 
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These results support two concepts. First, for a fixed grid density, as the 
order of the specified filter increases (and hence the number of specified frequency 
samples), the likelihood of the occurrence of a degeneracy increases. Secondly, for a 
fixed filter order, as the density of the grid on which samples are placed increases, 
the likelihood of a degeneracy occurring decreases. For filters of order 17 x 17 
or less, the probability of a degenerate case is minimal provided frequencies are 


specified to at least the nearest 0.05 radians. 


TABLE 3-1 
NO. OF DEGENERACIES OCCURRING IN 100 TEST CASES 







GRID DENSITY (RADIANS) 


FILTER ORDER RO 0.2 0.1 0.05 


+) ] 0 
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It is significant to note that for the above tests, frequency samples were cho- 
sen in an arbitrary manner. However, when designing a filter to approximate 
a certain desirable frequency response characteristic, the selection of points is 
likely to be done in an unconstrained, but not arbitrary, manner. Certain pat- 
terns in sample selection, such as selecting octant symmetric sample points where 
(w,W2); = (w2,u,);, have a greater tendency to introduce linear dependencies 
within the set of characteristic vectors than when sample points are selected at 


random. This concept is difficult to show in a controlled experiment; however, in 
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the course of this research, degeneracies arose when attempting filter designs much 


more often than the arbitrarily placed sample test results would indicate. 


Such a degeneracy did arise in the following design example. In this case, a 
lowpass filter is desired with a circularly symmetric frequency response such that 
it possesses unity gain in the passband region {0 < \/w? +w? < 0.27} and zero 
gain in the stopband region {0.47 < ye < a}. The region of support for 
the impulse response is specified to be 17 x 17, implying Af, = Af, = 8. Therefore, 
a total of R = (M, +1)(M, +1) = 81 frequency samples are to be specified. Figure 
3.3 shows the location of 81 points at which the desired frequency response is to 


be sampled initially. 
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Figure 3.3. Arrangement of Samples for Lowpass Filter 
Design, Case 1. 
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Note that the arbitrary-sampling method does not require any samples to be 
placed in the transition (don’t care) region. This allows the concentration of more 
samples in the passband and stopband regions. This particular choice of sample 
locations, however, was found to result in a degeneracy when a computer solution 
was attempted. 

Following this degenerate result, the sampling geometry was altered by slight 
displacement of the locations of six of the frequency samples. The new sampling 
geometry is illustrated in Figure 3.4, with the circled points indicating the sample 
locations which differ from the degenerate case. The resulting filter frequency 


response is shown in Figure 3.5. 
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Figure 3.4. Arrangement of Samples for Lowpass Filter 
Design, Case 2. 
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Figure 3.5. Frequency Response of Lowpass Filter 
Design, Case 2. 


Although the result may seem to be perfectly acceptable, the system of equa- 
tions involved in the solution can be considered to be ill-conditioned. Slight per- 
turbation of the position of a single sample point, as shown by the circled point in 
Figure 3.6, results in the unacceptable, excessively rippled response indicated by 
Figure 3.7. 

The above example demonstrates a negative aspect of the arbitrary-sampling 
approach. Although degeneracies, as defined, may seldom occur, ill-conditioned 
cases in which the determinant of matrix B in Eq. (3.12) is very small , but non- 
zero, may arise much more frequently. While specifying the frequencies w, and wy» 
at each sample point with a high degree of precision (analogous to placing samples 
on a very dense grid) may assure a very minimal probability of a degeneracy, it 


does not guarantee that the resulting system will not be ill-conditioned. 
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Figure 3.6. Arrangement of Samples for Lowpass Filter 
Design, Case 3. 


C. COMPUTATIONS 

A further negative aspect of this arbitrary-sampling method is the number 
of computations involved in determining filter parameters. While the uniform- 
sampling method of Chapter JI was shown to be reasonably efficient in computation, 
the arbitrary-sampling method involves the solution of a rather large system of 
linear equations. This solution requires the inversion of an R x R matrix. This 
inversion dominates the computation process. It is well known that inverting an 
NxWN matrix requires on the order of N° floating point operations, where a floating 
point operation consists of one multiplication and one addition. Since, for a filter 


with an N, x Nz impulse response region of support, 


N, -1 


(Sarva aah 


-(N, “FUNG Sea) 





R=(M, +1)(M +1) +1) 


er 
2 (3.18) 
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Figure 3.7. Frequency Response of Lowpass Filter 
Design, Case 3. 


the order of the number of floating point operations involved in the arbitrary- 


sampling design process can be represented as 
1 3 3 
| Cae ae Nes (3.19) 
For the case in which N, = Nz = N, this reduces to the following: 
Cae | Ne. (3.20) 
arb 64 


Thus, the arbitrary sampling method represents an approach which, while 
quite flexible, is computationally intensive and fails to provide guarantee of a 
unique solution. In the next chapter, an alternative nonuniform sampling approach 
is presented which, by placing some mild constraints on frequency sample location, 
always provides a unique design solution and is much more computationally effi- 


cient. 
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IV. A NOVEL 2-D NONUNIFORM FREQUENCY-SAMPLING APPROACH 





Chapters II and II] presented two fundamentally different approaches to 2—D 
FIR frequency-sampling filter design. The uniform-sampling method possessed no 
flexibility in the location of frequency samples, yet was computationally efficient in 
determining filter parameters. The arbitrary sampling method, however, had few 
constraints on sample location but was computationally intensive. The approach 
described in this chapter represents a compromise. It is more constrained than the 
arbitrary sampling method, but has much more freedom in sample placement than 
for the uniform-sampling method. While not as efficient as the IDF T approach, 
this method requires significantly fewer computations when compared with the 
arbitrary sampling method. Additionally, existence and uniqueness of a solution 
for the impulse response is guaranteed, provided simple constraints on sample 


locations are met. 


A. APPROACH 

For convenience, the design method introduced in this chapter will be referred 
to as the “nonuniform sampling” method, although the arbitrary sampling method 
uses nonuniformly spaced samples as well. This method borrows ideas from the 
other two approaches. As in the arbitrary-sampling case, quadrant symmetry is 
specified. As was shown in Chapter III, the frequency response of a 2—D zero-phase 


FIR filter with an N, x Nz impulse response region of support can be represented 


as 
Mee MM; 


Eh 5 | = » Se a(n, ,N2) COS W, 7, COS W2N2, (4.1) 


nya =—U ns —O 
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where 
Wigs >) ae 
= 2h(n1, 2) Ny 22 QO, MN = 0 
nea) 2hiny, nr 10 
4h(n,, nz) Ny oa QO, N2 - 0 


with M, = (N, — 1)/2 and M, = (N2 — 1)/2. In the arbitrary sampling method, 


(4.2) 


the sequence a(n,,n2) was represented as a column vector, the basis functions 
w;(w,,w2) were specified, and a system of equations of the form of Eq. (4.1) was 
solved directly. In the nonuniform sampling method, constraints placed on the 
sampling geometry result in Eq. (4.1) becoming a separable system. The system 
is then solved by a two-step process analogous to the row-column decomposition 
used in the uniform sampling method of Chapter II. The resulting systems of linear 
equations contain one-dimensional basis functions only; hence, linear dependencies 
associated with two-dimensional basis functions are not of concern. This results in 
assurance of a unique solution to the filter design problem. 

As in the arbitrary sampling method, samples are to be placed in region K = 
{(wi,W2) :0< w, < 7;0 < we < a} of the two-dimensional frequency plane. 
As mentioned above, certain constraints on the location of frequency samples are 
required. Suppose any (M, +1) distinct values of w, are chosen on the interval 
O<w, <a. These values are denoted as {w,, :k =0,1,...,M,}. Then for each 
W1,, Suppose any (M, +1) distinct values of w, are chosen such that 0 < w, <7. 
These values are denoted as {w;(k) :/ =0,1,...,M2}. Then a unique solution for 


a(n,,n,) can be obtained from 


(wi ~ Wei (k 2 5g a(n, Nz) COS Wy, 1 COS W2;(k) 9; (4.3) 
nmy=O nts—0 ° 
O<k<M,, 0<I1<M. 
Such a sampling geometry is shown in Figure 4.1 for the case M, = M, = 2. 


Notice that the selection of sample points is quite flexible. The values wo, 4,1, 
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and w,2 can be any arbitrary, distinct samples of w, in the region. Additionally, 
Wo0(k), w2,(k), and w22(k) represent any arbitrary, distinct samples of w2 for each 
k, k = 0,1,2. This is far less constrained than placing the samples on a uniform 


Cartesian grid. 


Ow, 
: @ ~22(0) 
@ “22(2) 
@ 21(0) 
@ w21(2) 
@ ~>2(1) 
@ ~2:(1) 
@ wo(1) 
e w20(0) 





Figure 4.1. Nonuniform Sampling Example. 


Given the arrangement of samples described above, the sequence a(n,,n2) can 


be found as follows. First, Eq. (4.1) is arranged as 


Mis aah 
H(w,,w2) = SS by a(n, m2) coswin; | cos wno. (4.4) 
ig Oi) =v 


A new function g(w,,n-) is defined as 


My, 


MO Oy = Ss a(n,,n2) cosw,n,. (4.5) 


n,=0 
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The frequency response can then be expressed as 


M 3 


HG) — = g(w1,N2) COS W2N2. (4.6) 


m4=—0 
For a particular value of w,, namely w,,, the response can be written as 


M 3 


(sone) — SS g(Wy_, Nz) COS WN. (4.7) 


m4=0 
Now, for each w2;(k), 1! =0,1,..., M2, the above expression becomes 


M 3 


H (wi, Wai(k)) = ww g(Wip, M2) COS Wo, {k) no. (4.8) 


a 0 
Thus, the first step in the row-column decomposition is to solve a system of (M, +1) 


linear equations of the form of Eq. (4.8). In matrix form, this becomes 


H (1%, W20 (k)) 1 coswoo({k)  ... cos Maw {k) g{w;,0) 

FH (w1~,We1(k)) 1 cosw2,(k) ... cos Mowe, (k) g(wy,, 1) 

H (14, Wem, (k)) 1 coswoy,(k) ... cosMawoy,(k) 4d Lg(wir,Me) . 
(4.9 


In vector-matrix notation, Eq. (4.9) is expressed as 
H, — Ag, 5 (4.10) 


where H, and g, are (M, +1) x1 column vectors, and A is an (M, +1) x (M, +1) 


square matrix. The vector g, is then computed using the relation 
g, =~ He (4,105) 


This entire process is performed for each k, k = 0,1,...,M,. Therefore, this step 
involves a total of (M, + 1) systems of (M2 + 1) linear equations each. 
Once the column vectors g, have been computed for all values of k, a matrix 


G of size (M2 + 1) x (M, +1) can be formed such that 
G= | Zo Si eee Em, ie (4.12) 
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Now the second step of the row-column decomposition can be performed. From 
Eq. (4.5), for a particular w4,, 


My, 


MOA) = yy a(n,,n2) cosw,, 7. (4.13) 


111 —0 
For a particular n,, a system of linear equations can be formed by evaluating Eq. 


(4.13) for each k, kK =0,1,...,M,. In matrix form, this becomes 


g(Wi0, 72) 1 COS Wig eee cos M, wo a(0, N2 ) 
g(wi1,%2) 1 COS W)1]1 eee cos Mw, a(1, 72) 
=]. , ; (4.14) 
G(Wim, >To) 1 coswijy, --- cosM,u,y,, a(M,, nz) 
Equation (4.14) can be expressed in vector-matrix notation as 
r,,, = Ba, (4.15) 


where I, is a (M, +1) x 1 column vector formed from the elements of the nth 
row of G, B is a (M, +1) x (M, +1) square matrix and a is a (M, +1) x 1 column 


vector. Thus a is found from the relation 
a — De lee (4.16) 


The vector a represents the values of a(n, ,n~) for a specified value of n.. To obtain 
the entire (M, + 1) x (M, + 1) sequence a(n,,n2), the process described by Eq. 
(4.16) is repeated for each value of n2, n. = 0,1,...,M.. Therefore, a total of 
(M, + 1) systems containing (M, + 1) linear equations each are involved in the 
second step of the method. Finally, the impulse response of the resulting filter can 
be evaluated using the relationship of Eq. (4.2). 

The following design example illustrates how this method can be employed. 


For this example, a 2—D circularly symmetric lowpass filter characteristic is desired 
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with unity gain in the passband region {(w,,w2) :.0 < V/w? +w? < 0.27} and zero 
gain in the stopband region {(w,,w2) : 0.47 < y oo < a}. The region of 
support for the impulse response is specified to be 17 x 17. The arrangement of 
frequency samples shown in Figure 4.2 is to be used. Note that the sampling is 
not uniform. When samples are placed in the transition region, a linear slope 
characteristic between the two bands is used in specifying the desired frequency 
response value at those samples, as was done in the uniform sampling approach. 
Contour and perspective plots of the frequency response for the resulting filter are 


shown in Figure 4.3. 
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Figure 4.2. Arrangement of Samples for 17 x 17 Lowpass Filter. 


One rather obvious extension of the nonuniform sampling method involves 
again specifying an N, x N. impulse response region of support, but in this case 
performing the interpolation in the first step of the process on w,. Suppose any 
(M, +1) distinct values of w. are selected on the interval 0 < w, < a. These are 


denoted as {w2,, ! = 0,1,...,M2}. For each value of w2,, any (M, + 1) distinct 
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Figure 4.3. Frequency Response of 17 x 17 Lowpass Filter. 
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values of w, are chosen on the interval 0 < w, < 7. These values are denoted as 
{wix(l), k= 0,1,...,M,}. Such an arrangement of selected samples is illustrated 
in Figure 4.4 for the case Af, = M, = 2. Now suppose H(w,,wz2) is specified at 


each selected sample (w,,(l),w2:) for k = 0,1,...,M, and! =0,1,...,M. 





Figure 4.4. Alternate Nonuniform Sampling Approach. 


Then there is a unique solution for a(n, ,n.) which can be obtained from 


M, 
H (Gye) oon) =» ye a(n,,n2)cos uw, (l)n, cos wn; 


m,y=On,=0 


Of ks M,, Usa a 


(4.17) 


The approach to the solution of the above problem and proof of its existence and 
uniqueness can be performed in a straightforward manner analogous to that of the 


basic approach. 
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It is worthy to note that the sampling geometry required by the uniform 
sampling method of Chapter II is, in reality, a special case of the general nonuni- 
form sampling form discussed here. Therefore, the results of the uniform sampling 
method can always be achieved by constraining the location of samples to the uni- 
form Cartesian grid. However, the method of solution is less efficient than the 


IDFT approach and so would probably not be used for this particular case. 


B. EXISTENCE AND UNIQUENESS 

A significant aspect of the nonuniform sampling method is that a unique solu- 
tion matching the specified frequency response at each sample point is guaranteed, 
provided samples are located as described. Proof of existence and uniqueness of a 
sequence h(n,,n2) which satisfies the interpolating conditions follows. 

The first step in the nonuniform sampling approach involves solution of (M, + 
1) systems of linear equations of the form of Eq. (4.9) or Eq. (4.10). A unique 
solution exists for each of these systems. Recall that Eq. (4.9) resulted from 
applying the interpolating conditions to Eq. (4.6). Through use of the Chebyshev 
polynomials, Eq. (4.6) can be expressed as 


M 3 


ial, ae, ye Cl, steomlicosieys) 35 (4.18) 


nz=0 
where a unique relation exists between the coefficients g(w,,n2) and c(w,,n). 


Through the subtitution z = cosw., the above becomes 


M3 
owen Clones cee (4.19) 
mz2=0 
The set {1,2,27,...,2%?} is known to be unisolvent on any interval [a,b] [Ref. 9: 


p. 31]. Therefore, a unique interpolation polynomial of the form of Eq. (4.19) can 


be found given any (M, + 1) distinct values of x. From the relation xz, = cos wy, 
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as long as distinct w2,(k) samples are specified on the interval {0 < w. < 7}, the 
corresponding z, will be distinct. Due to the unique relationship between Eas. 
(4.6), (4.18), and (4.19), the interpolation implied by Eq. (4.10) will yield a unique 
solution for g., kK = 0,1,...,M,. This implies that the matrix G defined in Eq. 


(4.12) will be unique. 


The second step of the nonuniform sampling approach involves solution of 
(M, +1) systems of linear equations of the form of Eq. (4.15). Each row of matrix 
B in Eq. (4.15) forms a vector of Chebyshev polynomials which form a unisolvent 
set. Therefore, B is invertible. Since a unique matrix G was found in the first step, 
the column vectors T,,, are unique for each nz = 0,1,..., M2. Therefore, a unique 
solution for a(n, ,n2) exists which satisfies the overall interpolation problem. From 
the relationship indicated by Eq. (4.2), this implies existence and uniqueness for 


the resulting impulse response h(n, ,nz2). 


[ll-conditioned results normally occur in interpolation problems involving one 
variable when one sample is located in close proximity to another. Since the 
nonuniform sampling method uses a two-step process involving one-dimensional 
interpolating functions, the same result can be expected to hold here. Provided 
samples are not placed too close together, the resulting filter design should not be 


ill-conditioned. 


C. COMPUTATIONS 


While the nonuniform sampling method retains less flexibility than the ar- 
bitrary sampling method of Chapter III, it is more computationally efficient for 
filters of the same order. As in the analysis for the arbitrary sampling case, in- 
verting an N x N matrix is assumed to require on the order of N° floating point 


operations and additions, and inversion of the matrix is assumed to dominate the 
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computations required in solving a system of linear equations. A brief analysis of 
the number of operations required for the nonuniform sampling method follows. 

For a filter with a specified impulse response region of support of size N, x N2, 
step one of the design process requires solution of a total of (M, + 1) systems of 
(M, + 1) linear equations each. The number of floating point operations involved 
in this step is therefore on the order of (M, + 1)(M, +1)°. 

Step two of the design process requires solution of a total of (M, +1) systems 
of (M, +1) linear equations each. The number of required floating point operations 
in this step, then, is on the order of (M, +1)(M, + 1)°. 

The total number of floating point operations for the entire nonuniform frequency- 
sampling 2-D FIR filter design method is the combined amount from the above 


two steps. The total number, therefore, is on the order of 
Cau = (Mi +1)(M2 +1) + (M +:1)°(M +1) 


= (A 41) (4) (ty) (tt) 
7 DD 2 2 2 


1 eam ; (4.20) 
= —(N, +1)(N. +1)? + —(N, +:1)°(N2 +1) 
16 16 
JL 
For the particular case in which N, = N, = N, Eq. (4.20) reduces to 
1 4 


Thus, the nonuniform sampling method is much more computationally effi- 
cient than the arbitrary sampling method, yet is still quite flexible while providing 


guarantee of a unique design solution. 


47 


V. A COMPARISON OF FREQUENCY-SAMPLING DESIGN METHODS 
THROUGH DESIGN EXAMPLES 


Up to this point the theory behind three fundamentally different approaches 
to the frequency-sampling design of 2—D FIR filters has been developed. With 
knowledge of these approaches, it is of interest to investigate the results of applying 
each of the methods to some practical design problems. In this chapter, some design 
problems are formulated and the results from each of the methods are compared 


in terms of performance, computations, and flexibility. 


A. APPROACH 

In order to compare the three basic design methods discussed in this thesis, 
three different design problems are described. To achieve a valid comparison, sev- 
eral parameters were fixed for each design. The region of support for the filter 
impulse response in all cases was specified to be 17 x 17. A square region of sup- 
port is most commonly used in 2—D filter design and avoids some of the problems 
discussed in the latter part of Chapter IV. The specified dimensionality of the im- 
pulse response is large enough to ensure reasonable results but not so large as to 
make implementation impractical. 

Each design method requires specification of a desired frequency response 
H,(w,,w,). For each design method, H,(w,,w2) is fixed, with unity gain through- 
out the passband, zero gain throughout the stopband, and a linear slope charac- 
teristic in the transition band(s). The frequency response value assigned to each 
sample point in the two-dimensional frequency domain is the value of Hz(w,,w2) 


at that point. 
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For each design problem, circularly symmetric filters are prescribed. Although 
more general symmetries are allowable for each method, circularly symmetric de- 


signs adequately illustrate the salient characteristics of the methods. 


A few other points are of note. For these design examples, optimimization 
techniques were not utilized in locating frequency samples for any of the methods. 
Frequency sample locations were selected strictly on a trial-and-error basis by the 
human designer for the arbitrary and nonuniform sampling methods. The resulting 
frequency responses of the filters are compared qualitatively in terms of passband 
and stopband ripple, transition band width, and circularity. No quantitative com- 
parisons are made in terms of best equi-ripple approximation to an ideal response 
characteristic, etc. It is well known that a better quality filter can usually be ob- 
tained by increasing the number of terms in the impulse response (i.e. increasing 
the filter order). In comparing frequency-sampling methods with the same support 
for the impulse response in each case, substantially better filters are not expected 
to result from one particular method; rather, it is of interest to investigate how 
design method flexibility can be used to enhance certain desirable characteristics 


in the resulting filter frequency response. 


B. DESIGN EXAMPLES 


Three different design problems are now presented. The three frequency- 
sampling design methods are applied to each and a comparative analysis of the 


resulting filter frequency responses is performed. 


1. Lowpass Filter No. 1 
For this example, a circularly symmetric 2—D lowpass filter frequency re- 
sponse is desired such that it has unity gain in the passband region {0 < \/w? +w? < 


0.47} and zero gain in the stopband region {0.6m < \/w? + w?2 < mr}. First,a 17x17 
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filter was designed to match the above characteristics using the uniform sampling 
method. Contour and perspective plots of the filter frequency response are shown 
in Figure 5.1. Next, a 17 x 17 filter was designed using the nonuniform sampling 
method with the arrangement of samples shown in Figure 5.2. The resulting fil- 
ter frequency response is illustrated in Figure 5.3. Finally, the arbitrary sampling 
method was applied to this design problem. The selection of samples was as shown 
in Figure 5.4, and the frequency response of the resulting filter can be seen in 
Figure 5.5. 

In comparing the results of the three methods, the quality of the result- 
ing filters varies slightly in each case but the arbitrary and nonuniform sampling 
designs appear no better than that of the uniform sampling case. The uniform 
sampling method required only 2 x 17° = 9826 floating point operations, while 
the nonuniform sampling method required on the order of 2 x 9* ~ 13,000 float- 
ing point operations. The arbitrary sampling method required on the order of 
1/64 x 17° ~ 380,000 floating point operations. In looking back at the desired 
characteristics of this filter, the passsband and stopband regions comprise fairly 
equal areas of the (w,,w2) plane. Sampling on a uniform Cartesian grid allows 
an adequate number of samples in all bands and hence little seems to be gained 
through nonuniform spacing of samples. 

2. Lowpass Filter No. 2 

For the second example, another circularly symmetric 2—D lowpass filter 
characteristic is desired. In this case, the filter is to have unity gain in the passband 
{0 < \/w? + w? < 0.27} and zero gain in the stopband {0.47 < \/w? +w? < zy}. 
This particular example was selected since it illustrates a case in which the desired 
passband encompasses a much smaller portion of the region K of the (w,,w2) plane 


than the stopband. 
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Figure 5.1. Frequency Response of Lowpass Filter Design No. 1, 
Uniform Sampling Method. 
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Figure 5.2. Arrangement of Samples for Lowpass Filter 
Design No. 1, Nonuniform Sampling Method. 
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Figure 5.3. Frequency Response of Lowpass Filter Design No. 1, 
Nonuniform Sampling Method. 
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Figure 5.4. Arrangement of Samples for Lowpass Filter 
Design No. 1, Arbitrary Sampling Method. 
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Figure 5.5. Frequency Response of Lowpass Filter Design No. 1, 


Arbitrary Sampling Method. 


The use of the uniform sampling method allows the placement of only four 
samples in the passband region. The resulting filter frequency response, shown in 
Figure 5.6, has a rather smooth passband and stopband characteristic, but the 
passband is smaller in area than desired and can be seen to have a rather squared 
shape in the contour plot. 

The specified selection of frequency samples for the nonuniform sampling 
method is illustrated in Figure 5.7. Note that five samples are now placed within 
the desired passband region. Figure 5.8 represents the frequency response of the 
resulting filter. In this case, the size of the resulting passband more closely matches 
the desired characteristic and the passband is clearly more circular in shape. These 
characteristics are gained, however, at the expense of introducing a greater amount 
of ripple in the stopband. 

Figure 5.9 shows a chosen arrangement of samples for the arbitrary sam- 
pling method. In this instance, six samples are now contained in the specified 
passband region while no samples are placed in the transition region. The result- 
ing filter frequency response is illustrated in Figure 5.10. 

Note that for the latter result, the passband is fairly flat and quite circular, 
and the roll-off is quite sharp. This occurs since the flexibility of the design allows 
frequency samples to be placed along the perimeters of the specified passband 
and stopband regions for the latter two design methods. Because of this, the two 
methods which allow flexibility in sampling can be used to provide greater control 
over the resultant width of the transition band. Of course, this is achieved at the 
expense of increased ripple. Since the arbitrary and nonuniform sampling methods 
do not require sampling on a Cartesian grid, they are generally better able to 


conform to non-rectangularly shaped contours such as circles. 
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Figure 5.6. Frequency Response of Lowpass Filter Design No. 2, 
Uniform Sampling Method. 
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Figure 5.7. Arrangement of Samples for Lowpass Filter 


Design No. 2, Nonuniform Sampling Method. 
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Figure 5.8. Frequency Response of Lowpass Filter Design No. 2 
Nonuniform Sampling Method. 
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Figure 5.9. Arrangement of Samples for Lowpass Filter 
Design No. 2, Arbitrary Sampling Method. 


60 


Hi@ tain) 





Figure 5.10. Frequency Response of Lowpass Filter Design No. 2, 
Arbitrary Sampling Method. 
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3. Bandpass Filter 

The final design example illustrates how each of the frequency-sampling 
techniques can be used in the design of a 2—D circularly symmetric bandpass filter. 
Here, the desired filter has unity gain in the passband region {0.47 < Jw? + w? =< 
0.67}, and zero gain in the stopbands defined by {0 < \/w? +w? < 0.27} and 
Orda << Jw? + w? < a}. The filter resulting from application of the uniform 
sampling method has the frequency response illustrated in Figure 5.11. For the 
nonuniform sampling method, the arrangement of samples shown in Figure 5.12 
was used. The resulting filter frequency response is shown in Figure 5.13. Figure 
5.14 illustrates the arrangement of samples selected for the arbitrary sampling 
method. The resulting filter has the frequency response of Figure 5.15. 

In comparing these results, the nonuniform sampling method produces a 
somewhat flatter and wider passband than the uniform sampling method, but re- 
sults in some notable ripples in the outer stopband which are not found in the 
uniform sampling case. The passband for the arbitrary sampling design is even 
more wide and flat, but at least for this particular selection of samples, has unac- 
ceptably high ripple in both stopband regions. 

The point to be made in comparing these methods is not that the nonuni- 
form or arbitrary sampling methods necessarily produce better filters, but rather 
that they possess more flexibility and hence can be used to enhance a particular 
desirable feature at the expense of another. The advantage gained through design 
flexibility may well outweigh the cost of additional computations for many design 


problems. 
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Figure 5.11. Frequency Response of Bandpass Filter Design, 


Uniform Sampling Method. 
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Figure 5.12. Arrangement of Samples for Bandpass Filter 
Design, Nonuniform Sampling Method. 
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13. Frequency Response of Bandpass Filter Design, 
Nonuniform Sampling Method. 
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Figure 5.14. Arrangement of Samples for Bandpass Filter 
Design, Arbitrary Sampling Method. 
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Figure 5.15. Frequency Response of Bandpass Filter Design, 


Method. 


Arbitrary Sampling 
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VI. CONCLUSIONS 


In this thesis, various approaches to the frequency-sampling design of two- 
dimensional FIR digital filters were investigated. The traditional frequency-sampling 
approach, referred to as the uniform sampling method here, involves specifying the 
frequency response at sample points located at the vertices of a uniform Cartesian 
grid in the (w,,w,) plane. Through use of the IDFT, the method is quite compu- 
tationally efficient, and can be made even more efficient if an FFT is used. The 
method possesses no flexibility in matching design parameters, however. Another 
frequency-sampling approach arises from specifying the frequency response at ar- 
bitrarily selected distinct samples in the (w,,w2) plane and then solving a large 
system of linear equations. This method, referred to as the arbitrary sampling 
method in this thesis, possesses great flexibility but is computationally intensive 
and prone to the occurrence of degeneracies in which a unique solution does not ex- 
ist. Finally, a new frequency-sampling approach, termed the nonuniform sampling 
method, was introduced. This method somewhat constrains the location of fre- 
quency samples, but still provides some free parameters and is reasonably efficient 
in computation. 

The results of this thesis reflect how computational efficiency and design flexi- 
bilty relate to the number of constraints on sample location. The methods described 
here which impose more constraints are more computationally efficient. Methods 
with less constraints were found to be more flexible in matching a desired frequency 
response. 

A problem posed by both the arbitrary and nonuniform sampling methods is 


how to properly select the samples in order to achieve an acceptable filter. Various 
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algorithms for the design of 2-D optimal filters have been developed which involve 
the solution of linear equations, analogous to those for the arbitrary sampling 
method, in an iterative approximation process. Although they have been shown to 
converge, the algorithms necessarily take special measures to deal with degeneracies 
and are very computationally intensive. [Refs. 9, 10] 

The nonuniform sampling method warrants a closer look because of its rela- 
tive computational efficiency. Due to constraints on the sampling geometry, the 
nonuniform sampling method is probably not well-suited to the design of globally 
optimal filters. If an efficient algorithm can be developed which will locate a set of 
points satisfying the non-uniform method constraints, and if this set guarantees a 
filter which, while not necessarily optimal, nevertheless is of acceptable quality, the 
non-uniform sampling method developed in this thesis may prove to be of great 


significance. 
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APPENDIX 


A 2-D FIR FILTER DESIGN BASED ON A BIVARIATE EXTENSION 
OF NEWTON’S POLYNOMIAL INTERPOLATION FORMULA 


In the course of this thesis research, an extension of Newton’s polynomial 
interpolation method in two variables was studied and applied to the design of two- 
dimensional digital filters. Such an extension of Newton’s method was proposed by 
Diamessis [Ref. 11]. Its application as a 2-D FIR filter frequency-sampling design 
method was presented at the ICASSP 87 conference [Ref. 12] and is summarized 


here. 


A. 2-D NEWTON INTERPOLATION 
The bivariate Newton interpolation problem for the 2-D quadratic case can 
be described as follows. Given 


1) a set of points in the (z,y) plane: 


S, = {(20,40), (20.91), (Zo, Yo), (21, Yo), (21,91), (22 Yo) } (A 


2) a set of values at those points: 


SS f (Zo, yo) = foo, f (Zorn) — ie Pa su, = 
FEL Hen) eho flti,m) 2 fii, f(z2,yo) = fro Ge 


3) a quadratic 2-D polynomial: 
p(z,y) = Apo + Ao Ft Go y+ Agr +41, FY + Ap2y” (A 7 3) 
find 


A Tr 
a= 290 > 410; 401» 220, 411; 42 | (A — 4) 
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such that the following interpolating conditions are satisfied 


P(Zo,Yo) = foo. PlZo.y1) = for, PlZo0,¥2) = foe 7 eer) 
P(Z1,Yo) = fio, P(Zi.wi)=fir, Plr2,¥o) = foo 
If the polynomial of Eq. (A-3) is rewritten in the form 
Pie i—teyy Cro (2 =o) fh Cony = Uo) 
+ Coo (x — Zo)(z — 2) + C11 (2 — Zo) (y — Yo) (A — 6) 
+corly—y)l(y—m), 
then the interpolation problem involves finding the coefficient vector 
Geee lenanc cane une ACT Cae | (A — 7) 


instead of a. The system of equations formed by application of Eq. (A-2) to Eq. 


(A-6) is 
1 0 @) 0 0 O Coo 
1 (x, = Zo) 0 0 0 O Cio 
1 0 (y1 — yo) 0 0 0 col 
1 (z2- Zo) 0 (zz ~ Zo) (z2 - 2) 0 0 C20 
1 (z1—2Z0) (y1 — yo) 0 (z1 — Zo)(y1 — yo) 0 C11 
1 0 (y2 — yo) 0 0 (y2 — yo)(ye -—yi)d Leo 
foo 
fio 
_ | for 
f20 
via 
fo2 
(A —8) 
or, In vector-matrix notation, 
Nop c—T. (A 5 9) 


Since N.p is lower triangular, c can be solved for recursively. Also, Np is 
non-singular provided the z; are distinct and the y; are distinct. Hence a unique 
solution for c can always be found. The coefficient vector a is related to c by the 
relation 


c= U.pa, (A — 10} 


where U,p is an upper triangular matrix which transforms the Newton coefficient 
vector c into the coefficient vector a. Therefore, a can be found through solution 


of the two triangular systems indicated by Eqs. (A-9) and (A-10). 


B. NEWTON FILTER DESIGN METHOD 

The selection of frequency samples for the Newton design method follows from 
the two variable Newton interpolation method discussed above. This approach is 
a 2-D extension of a one-dimensional filter design proposed by Schuessler |Ref. 13: 
p. 257]. If a square region of support on the first quadrant of extent N x N is 
initially specified for the filter impulse response, the frequency response of such a 
2—-D FIR filter with linear phase can be expressed as 


N-1N-1 


H(w,, We) oy \ h(nians eee: cm aac (A — 11) 


my=On,=0 
A total of (M+1)(M + 2)/2 frequency samples are then chosen on region K of the 
(w,,W2) plane, where M = (N —1)/2 and K = {(w,,w2),0 < wy, < 730 < we < TH. 


The selected samples are of the form 
{((wr,, Wo, ), Ke = 0,1,.2.;4kp) ee On eee (A — 12) 


A support set of the above form is referred to as triangular. Note that the samples 
need not be uniformly spaced. Since there is no constraint such that w,;, <w,;, < 


w;, °+*, the sampling geometry is not necessarily of a triangular shape on region 
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kK. Two such support sets are shown in Figure A.1. Given a specified value for 
H(w,,,W2,) at each sample location, the above support set ensures existence of a 
unique solution to the design problem. 


The Newton design imposes the quadrant symmetry condition 
wart) =n — 1 n-ne ae — 1 — 1). (A — 13) 


Application of Eq. (A-13) to Eq. (A-11) results in the expression 


M-1 
Hp (wi, w2) =h(M,M)+2 S— h(M,nz)Tu —n, (cos w2) 


ni (A — 14) 


M-1 2M 


ey, SS S_ A (n1,N2)Tu —n, (cos wi)Ty —n, (coswe), 


ny=On,=—0 
where T,, (xz) is the nth degree Chebyshev polynomial in z and Hy(w,,w-) is the 


zero-phase frequency response related to H(w,,w,) by 
id (ne, | =e 7 emt ee ty (ae). (A — 15) 
Now, if the impulse response region of support is restricted such that 
h(ni,n2)=0, |n, -M|+ |n.-—M|>M, (A — 16) 
and if the transformations 
cosw, =1—2zr and cosw,=1-2y (A — 17) 


are applied, then the zero-phase frequency response of Eq. (A-14) is mapped into 


a bivariate polynomial of the Newton form 


+—1 ae 


H(z, y) eo> co LT (cx — 2x, + Dew Tt (y— y,) 


il 7=1 


(A — 18) 
Nie t—j-1 he 

+) GH as Il @=5,) Ii@=a} 
= p=0 2 
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Lo Ly L5 og 


Yo 


Y 


Zo D2 Ly t 
(b) 
Figure A.1. Possible Sample Arrangements for Bivariate Newton 


Interpolation, Quadratic Case. 
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The design problem involves first mapping the specified frequency samples 
(w,, ,W2, ) into samples (z,,y,) through the relations of Eq. (A-17), and then using 
the Newton interpolation method previously described to find the coefficients c¢;,. 
By applying Eq. (A-17) to Eq. (A-18), and through the substitution 


COS W; = —_— «= > a 


a filter structure of the form 


M to 1 
Bee es 1 Z = 

H(z, , 22) =Coo 2, ™ z, 2 2, 60%; - ee = Il-70 — Pipé siete i *) 

C= p= 

M j-1 1 
~ (he) - —2 

fc) ines []- Gt — Monza! +25 ) 
j= 1 q= 


M t-1 


1—2 7 = 1 
t-j7-1 1 3-1 1 

[] -70- mer! +277) [] -¢0 -meaez' +23")}, 
p=O0 q=0 


(A — 19) 


results, where 4,, = 2cosw,, and U2, = 2coswa,. 


C. DISCUSSION 

The advantages hoped to be gained by design of the Newton filter were com- 
putational efficiency through the recursive nature of the solution and use of the 
permanence property {Ref. 11]. The latter property allows the specification of 
additional sample values for an interpolation problem without having to redeter- 
mine previously computed coefficients. These advantages, however, apply only if 
the filter is implemented using the structure implied by Eq. (A-19), where the co- 
efficients c;, represent the gains. Such a structure, however, requires more delays 
than a direct form structure of the same order which uses the impulse response 


coefficients h(n,,n~) for the filter gains, and hence is not very attractive. 
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An alternate procedure might be to use the Newton method to find the im- 
pulse response coefficients and implement the filter in direct form. Obtaining the 
impulse response coefficients of the filter from the set of coefficients c;; requires 
solution of two triangular system analogous to those of Eqs. (A-9) and (A-10). 
Carrying out this two-step process, solving first for c and then for h(n,,n2) to 
get a direct form implementation, is equivalent to solving a system of equations 
via LU-decomposition. This particular system is the system of linear equations 
resulting from the direct application of the interpolating conditions to Eq. (A-11), 
where the limits of the summation in Eq. (A-11) are altered to reflect the Newton 
filter region of support. This is the same approach used in the arbitrary sam- 
pling method of Chapter III, if the matrix inversion is performed through use of 
LU-decomposition, and therefore cannot be considered any more computationally 


efficient. 


The Newton method may be less significant in terms of computational effi- 
ciency or practical filter implementation than it is for studying the formulation 
of the two-dimensional frequency-sampling design problem. In Chapter III, it was 
shown that existence and uniqueness of a solution to the interpolation problem 
for arbitrarily placed samples and a fixed impulse response region of support was 
not guaranteed for the two-dimensional case. With specification of a triangular 
support set and a polynomial of the form of Eq. (A-6), the Newton method pro- 
vided a unique solution. If another, less constrained choice of samples is specified 
and if the form of the polynomial is altered to include basis functions such that 
the interpolating problem becomes a lower triangular system of equations akin to 
Eq. (A-8) with a non-zero main diagonal, then a unique interpolating polynomial 
of the new form will exist. Since the specified basis functions for the polynomial 


interpolation problem directly relate to the region of support for the impulse re- 
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sponse in the filter design problem, the Newton method can be used to illustrate 
a connection between the particular sampling geometry and an impulse response 


region of support which will ensure a unique solution. 


(es 


10: 


6) 


12: 
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